Supplement 3:
Supplementary information for McGowan et al. “Linking monitoring and data analysis to predictions and decisions for a range-wide eastern Black Rail status assessment,” including simulation model code for the occupancy simulation modeling to predict future status.

### Occupancy projection model for Eastern Black Rail SSA
###  Code drafted by C McGowan initiated on 1/3/2018, last edited 10/18/2019

yrs = 87 #number of years simulated
reps=5000  # number of replicates

####South East CP
Nocc = matrix(0,reps,yrs) #object to store simualtion output
Pocc=matrix(rbeta(reps,9,91),reps,1)
gPpersi = matrix(rbeta(reps,99,1),reps,1) # object to store replicate mean persistence prob
gSDppi = matrix(gPpersi*.1,reps,1) # object to store replicate persistence prob SD, SD are calculated here using a 10% coffecient of Variation
iPpersi = matrix(rbeta(reps,51,49),reps,1) # object to store replicate mean persistence prob
iSDppi = matrix(iPpersi*.1,reps,1) # object to store replicate persistence prob SD, SD are calculated here using a 10% coffecient of Variation
bPpersi = matrix(rbeta(reps,43,57),reps,1) # object to store replicate mean persistence prob
bSDppi = matrix(bPpersi*.1,reps,1) # object to store replicate persistence prob SD, SD are calculated here using a 10% coffecient of Variation

Pperst = matrix(0,reps,yrs) #object to store persistence probabilities for each year and rep
q=matrix(0,reps,yrs)
Tocc = 5000
Ty=matrix(Tocc,reps,yrs)
Hl=matrix(0,reps,yrs)
hlr=0.12 #change this for habitat loss scenarios
c=0.03
col=matrix(rbeta(reps*yrs,100*c,100*(1-c)),reps,yrs)
Nocci=matrix(0,reps,1) #Randomize the intial number of sites occupied using a unifor distribtion between 25 and 200 
ph=0.012 #change this for habitat quality scenarios
phab=matrix(rbeta(reps*yrs,100*ph,100*(1-ph)),reps,yrs)
Por=matrix(0,reps,yrs)

###Mid-Atlantic
Nocc2 = matrix(0,reps,yrs) #object to store simualtion output
Pocc2=matrix(rbeta(reps,9,91),reps,1)
gPpersi2 = matrix(rbeta(reps,99,1),reps,1) # object to store replicate mean persistence prob
gSDppi2 = matrix(gPpersi2*.1,reps,1) # object to store replicate persistence prob SD, SD are calculated here using a 10% coffecient of Variation
iPpersi2 = matrix(rbeta(reps,51,49),reps,1) # object to store replicate mean persistence prob
iSDppi2 = matrix(iPpersi2*.1,reps,1) # object to store replicate persistence prob SD, SD are calculated here using a 10% coffecient of Variation
bPpersi2 = matrix(rbeta(reps,43,57),reps,1) # object to store replicate mean persistence prob
bSDppi2 = matrix(bPpersi2*.1,reps,1) # object to store replicate persistence prob SD, SD are calculated here using a 10% coffecient of Variation

Pperst2 = matrix(0,reps,yrs) #object to store persistence probabilities for each year and rep
q2=matrix(0,reps,yrs)
Tocc2 = 3500
Ty2=matrix(Tocc2,reps,yrs)
Hl2=matrix(0,reps,yrs)
hlr2=0.12
c2=0.03
col2=matrix(rbeta(reps*yrs,100*c2,100*(1-c2)),reps,yrs)
Nocci2=matrix(0,reps,1) #Randomize the intial number of sites occupied using a unifor distribtion between 25 and 200 
ph2=0.012
phab2=matrix(rbeta(reps*yrs,100*ph2,100*(1-ph2)),reps,yrs)
Por2=matrix(0,reps,yrs)

###South west
Nocc3 = matrix(0,reps,yrs) #object to store simualtion output
Pocc3=matrix(rbeta(reps,24,76),reps,1)
gPpersi3 = matrix(rbeta(reps,85,15),reps,1) # object to store replicate mean persistence prob
gSDppi3 = matrix(gPpersi3*.1,reps,1) # object to store replicate persistence prob SD, SD are calculated here using a 10% coffecient of Variation
iPpersi3 = matrix(rbeta(reps,66,34),reps,1) # object to store replicate mean persistence prob
iSDppi3 = matrix(iPpersi3*.1,reps,1) # object to store replicate persistence prob SD, SD are calculated here using a 10% coffecient of Variation
bPpersi3 = matrix(rbeta(reps,41,59),reps,1) # object to store replicate mean persistence prob
bSDppi3 = matrix(bPpersi3*.1,reps,1) # object to store replicate persistence prob SD, SD are calculated here using a 10% coffecient of Variation

Pperst3= matrix(0,reps,yrs) #object to store persistence probabilities for each year and rep
q3=matrix(0,reps,yrs)
Tocc3 = 3500
Ty3=matrix(Tocc3,reps,yrs)
Hl3=matrix(0,reps,yrs)
hlr3=0.12
c3=0.13
col3=matrix(rbeta(reps*yrs,100*c3,100*(1-c3)),reps,yrs)
Nocci3=matrix(0,reps,1) #Randomize the intial number of sites occupied using a unifor distribtion between 25 and 200 
ph3=0.0002
phab3=matrix(rbeta(reps*yrs,100*ph3,100*(1-ph3)),reps,yrs)
Por3=matrix(0,reps,yrs)

###Great Plains
Nocc4 = matrix(0,reps,yrs) #object to store simualtion output
Pocc4=matrix(rbeta(reps,12,88),reps,1)
gPpersi4 = matrix(rbeta(reps,73,27),reps,1) # object to store replicate mean persistence prob
gSDppi4 = matrix(gPpersi4*.1,reps,1) # object to store replicate persistence prob SD, SD are calculated here using a 10% coffecient of Variation
iPpersi4 = matrix(rbeta(reps,42,68),reps,1) # object to store replicate mean persistence prob
iSDppi4 = matrix(iPpersi4*.1,reps,1) # object to store replicate persistence prob SD, SD are calculated here using a 10% coffecient of Variation
bPpersi4 = matrix(rbeta(reps,23,77),reps,1) # object to store replicate mean persistence prob
bSDppi4 = matrix(bPpersi4*.1,reps,1) # object to store replicate persistence prob SD, SD are calculated here using a 10% coffecient of Variation

Pperst4 = matrix(0,reps,yrs) #object to store persistence probabilities for each year and rep
q4=matrix(0,reps,yrs)
Tocc4 = 3000
Ty4=matrix(Tocc4,reps,yrs)
Hl4=matrix(0,reps,yrs)
hlr4=0.11
c4=0.03
col4=matrix(rbeta(reps*yrs,100*c4,100*(1-c4)),reps,yrs)
Nocci4=matrix(0,reps,1) #Randomize the intial number of sites occupied using a unifor distribtion between 25 and 200 
ph4=0.012
phab4=matrix(rbeta(reps*yrs,100*ph4,100*(1-ph4)),reps,yrs)
Por4=matrix(0,reps,yrs)

for(i in 1:reps){  ##Loop to create replicates
  #SouthEast
  Nocci[i]=rbinom(1,Tocc,Pocc[i])
  #Mid-Atl
  Nocci2[i]=rbinom(1,Tocc2,Pocc2[i])
  #SouthWest
  Nocci3[i]=rbinom(1,Tocc3,Pocc3[i])
  #GP
  Nocci4[i]=rbinom(1,Tocc4,Pocc4[i])
  for(j in 1:yrs){  ##Loop to create time lines
    #South east  
    q[i,j]=round(runif(1,0.5,3.5)) 
    if(q[i,j]==1) Pperst[i,j]=rbeta(1,100*gPpersi[i],100*(1-gPpersi[i]))-phab[i,j] #Drawing a randomized value for persistence probability
    if(q[i,j]==2) Pperst[i,j]=rbeta(1,100*iPpersi[i],100*(1-iPpersi[i]))-phab[i,j]
    if (q[i,j]==3) Pperst[i,j]=rbeta(1,100*bPpersi[i],100*(1-bPpersi[i]))-phab[i,j]
    Hl[i,j]=rbeta(1,100*hlr,100*(1-hlr))
    if(j>1)Ty[i,j]=round(Ty[i,j-1]*(1-Hl[i,j-1]))
    Nocc[i,1]=Nocci[i]#setting initial number of sites occupied
    if(j>1)Nocc[i,j]=rbinom(1,Nocc[i,j-1],Pperst[i,j])+(rbinom(1,Ty[i,j],col[i,j])) #Markovian process for future/simulated number of sites occupied
    Por[i,j]=Nocc[i,j]/Nocci[i]
  
    #Mid-Atl  
    q2[i,j]=round(runif(1,0.5,3.5)) 
    if(q2[i,j]==1) Pperst2[i,j]=rbeta(1,100*gPpersi2[i],100*(1-gPpersi2[i]))-phab2[i,j] #Drawing a randomized value for persistence probability
    if(q2[i,j]==2) Pperst2[i,j]=rbeta(1,100*iPpersi2[i],100*(1-iPpersi2[i]))-phab2[i,j]
    if (q2[i,j]==3) Pperst2[i,j]=rbeta(1,100*bPpersi2[i],100*(1-bPpersi2[i]))-phab2[i,j]
    Hl2[i,j]=rbeta(1,100*hlr2,100*(1-hlr2))
    if(j>1)Ty2[i,j]=round(Ty2[i,j-1]*(1-Hl2[i,j-1]))
    Nocc2[i,1]=Nocci2[i]#setting initial number of sites occupied
    if(j>1)Nocc2[i,j]=rbinom(1,Nocc2[i,j-1],Pperst2[i,j])+(rbinom(1,Ty2[i,j],col2[i,j])) #Markovian process for future/simulated number of sites occupied
    Por2[i,j]=Nocc2[i,j]/Nocci2[i]
  
    #SouthWest
    q3[i,j]=round(runif(1,0.5,3.5)) 
    if(q3[i,j]==1) Pperst3[i,j]=rbeta(1,100*gPpersi3[i],100*(1-gPpersi3[i]))-phab3[i,j] #Drawing a randomized value for persistence probability
    if(q3[i,j]==2) Pperst3[i,j]=rbeta(1,100*iPpersi3[i],100*(1-iPpersi3[i]))-phab3[i,j]
    if (q3[i,j]==3) Pperst3[i,j]=rbeta(1,100*bPpersi3[i],100*(1-bPpersi3[i]))-phab3[i,j]
    Hl3[i,j]=rbeta(1,100*hlr3,100*(1-hlr3))
    if(j>1)Ty3[i,j]=round(Ty3[i,j-1]*(1-Hl3[i,j-1]))
    Nocc3[i,1]=Nocci3[i]#setting initial number of sites occupied
    if(j>1)Nocc3[i,j]=rbinom(1,Nocc3[i,j-1],Pperst3[i,j])+(rbinom(1,Ty3[i,j],col3[i,j])) #Markovian process for future/simulated number of sites occupied
    Por3[i,j]=Nocc3[i,j]/Nocci3[i]
 
    #GP
    q4[i,j]=round(runif(1,0.5,3.5)) 
    if(q4[i,j]==1) Pperst4[i,j]=rbeta(1,100*gPpersi4[i],100*(1-gPpersi4[i]))-phab4[i,j] #Drawing a randomized value for persistence probability
    if(q4[i,j]==2) Pperst4[i,j]=rbeta(1,100*iPpersi4[i],100*(1-iPpersi3[i]))-phab4[i,j]
    if (q4[i,j]==3) Pperst4[i,j]=rbeta(1,100*bPpersi4[i],100*(1-bPpersi4[i]))-phab4[i,j]
    Hl4[i,j]=rbeta(1,100*hlr4,100*(1-hlr4))
    if(j>1)Ty4[i,j]=round(Ty4[i,j-1]*(1-Hl4[i,j-1]))
    Nocc4[i,1]=Nocci4[i]#setting initial number of sites occupied
    if(j>1)Nocc4[i,j]=rbinom(1,Nocc4[i,j-1],Pperst4[i,j])+(rbinom(1,Ty4[i,j],col4[i,j])) #Markovian process for future/simulated number of sites occupied
    Por4[i,j]=Nocc4[i,j]/Nocci4[i]
  }}

pdf("scenario1.pdf")
#SouthEast
Mn1=apply(Nocc,2,median)
Mp1=apply(Por,2,median)
LbMp1=apply(Por,2,quantile, probs = c(0.025))
UbMp1=apply(Por,2,quantile, probs = c(0.975))

plot(Mp1, ylab="Proportion of sites remaining occupied", xlab="year",ylim=c(0,2),xlim=c(0,yrs))
lines(Mp1,lty=1) 
lines(LbMp1, lty=2)#dashed line for median sites occupied for model 1
lines(UbMp1, lty=2)
title("Scenario 1, South East Coastal Plain")
legend("topright",c("Median", "C.I."), lty=c(1,2))

#Mid-Atl
Mn2=apply(Nocc2,2,median)
Mp2=apply(Por2,2,median)
LbMp2=apply(Por2,2,quantile, probs = c(0.025))
UbMp2=apply(Por2,2,quantile, probs = c(0.975))

plot(Mp2, ylab="Proportion of sites remaining occupied", xlab="year",ylim=c(0,2),xlim=c(0,yrs))
lines(Mp2,lty=1) 
lines(LbMp2, lty=2)#dashed line for median sites occupied for model 1
lines(UbMp2, lty=2)
title("Scenario 1, Mid-Atlantic Coastal Plain")
legend("topright",c("Median", "C.I."), lty=c(1,2))
#South-west
Mn3=apply(Nocc3,2,median)
Mp3=apply(Por3,2,median)
LbMp3=apply(Por3,2,quantile, probs = c(0.025))
UbMp3=apply(Por3,2,quantile, probs = c(0.975))

plot(Mp3, ylab="Proportion of sites remaining occupied", xlab="year",ylim=c(0,2),xlim=c(0,yrs))
lines(Mp3,lty=1) 
lines(LbMp3, lty=2)#dashed line for median sites occupied for model 1
lines(UbMp3, lty=2)
title("Scenario 1, Southwest Coastal Plain")
legend("topright",c("Median", "C.I."), lty=c(1,2))
#GP
Mn4=apply(Nocc4,2,median)
Mp4=apply(Por4,2,median)
LbMp4=apply(Por4,2,quantile, probs = c(0.025))
UbMp4=apply(Por4,2,quantile, probs = c(0.975))

plot(Mp4, ylab="Proportion of sites remaining occupied", xlab="year",ylim=c(0,2),xlim=c(0,yrs))
lines(Mp4,lty=1) 
lines(LbMp4, lty=2)#dashed line for median sites occupied for model 1
lines(UbMp4, lty=2)
title("Scenario 1, Great Plains")
legend("topright",c("Median", "C.I."), lty=c(1,2))
Supplement to McGowan et al. (2020) – Endang Species Res XX:XXX-XXX – https://doi.org/10.3354/esr01063